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Colloidal systems observed in video microscopy are often analysed using the displacements cor¬ 
relation matrix of particle positions. In non-thermal systems, the inverse of this matrix can be 
interpreted as a pair-interaction potential between particles. If the system is thermally agitated, 
however, only an effective interaction is accessible from the correlation matrix. We show how this 
effective interaction differs from the non-thermal case by comparing with high-statistics numerical 
data from hard-sphere crystals. 


1 Introduction 

The analysis of the displacement-displacement correlation func¬ 
tion between many particles has become a standard tool used on 
experimental data 1-5 acquired by video and confocal microscopy 
as well as on numerical data . 6-9 It is much, but not exclusively, 
used by researchers interested in jammed and glassy dynamics of 
soft particles. In practice, one first determines the reference po¬ 
sitions of particles and then records the correlations of displace¬ 
ments with respect to these reference positions for all particles. 
The averages are done either over time, using particle trajectories, 
or over an ensemble, using Monte-Carlo sampling. The result is a 
large correlation matrix G between all displacements of all parti¬ 
cles, which is then inverted to give another matrix D. 

This matrix D is often interpreted as an effective interaction ma¬ 
trix which characterises oscillations of the particles around their 
reference positions. In that context, D is called the dynamical 
matrix . 7,10 This interpretation has its origins in the analysis of 
non-thermal systems where massive particles are connected by 
harmonic or other interaction potentials. The dynamical matrix 
is then proportional to the Hessian matrix H of second deriva¬ 
tives of the total interaction potential 7 (with the particle mass 
and with kT as proportionality factors, which will play no role in 
the following). When activated by a small external stimulus, the 
response of the system is a mixture of oscillations which are given 
by the eigenfrequencies and eigenvectors of the dynamical matrix. 
One condition for this analysis to remain true is that the oscilla¬ 
tion amplitudes are small, so that the dynamical matrix does not 
depend on the displacements. With the masses of the particles 
and the measured correlation matrix G as input, one can obtain 
information about the interaction potential between the particles. 
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The protocol to record data and to invert the matrix G works 
equally well on thermal as on non-thermal systems. In contrast 
to the Hessian matrix H , which contains the static potential in¬ 
teractions, the matrix G _1 contains also all the thermal contribu¬ 
tions, including entropic. How is the dynamical matrix changed 
by these contributions? To what extent can we still infer the in¬ 
teraction between particles? For example, if the true interaction 
acts only between nearest neighbours, the same is true for H , but 
does it also hold for the effective interactions obtained from G _1 ? 

Between these two limits of completely non-thermal H and 
fully thermal G _1 , a plethora of intermediate levels can be identi¬ 
fied, which account for more or fewer aspects of the thermal dis¬ 
tribution and which may serve for more or less reasonable models 
for the observed system. We refer to Henkes et al 7 for a review 
of some of these intermediate levels and the connections between 
them, such as the harmonic approximation or the shadow system. 
As we mention models for the system, it is important to note that 
the Hessian matrix with which we compare G -1 may derive from 
modelling potentials very different from the true ones. 

The aim of the present paper is to establish the difference be¬ 
tween the two extreme cases, which are the non-thermal case 
(Hessian matrix H ) and the fully thermal matrix G -1 obtained 
from observation data. We here want to see the largest difference 
and therefore choose as system of observation a crystal of hard 
spheres which collide elastically. It is usually said that thermal 
(i entropic ) effects are maximal in the hard-sphere system. We also 
have to choose a model system of which we calculate the Hessian 
matrix. We here take the most general form of pairwise, central 
potential interactions. 

Beside its interpretation as effective interactions, the ma¬ 
trix G -1 - and more directly the correlation matrix G - serve also 
to extract elastic moduli from thermal fluctuations. Because the 
elastic moduli are macroscopic quantities of the whole system, 
this measurement involves only large-wavelength properties of 
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the effective interactions contained in G _1 . Of course, the same 
procedure can be done in the non-thermal case using H instead of 
G _1 . One advantage of the elastic moduli is that their calculation 
is not limited to smooth pairwise potentials but applies to other 
types of interaction as well, such as strictly hard spheres. We ex¬ 
pect to see here another discrepancy between the thermal and 
the non-thermal cases: When one searches the long-wavelength 
limit of the Hessian matrix, one obtains the so-called Cauchy rela¬ 
tions which reduce the number of independent elastic constants 
in materials below the number expected from symmetry consider¬ 
ations. For example, in two-dimensional hexagonal crystals and 
in three-dimensional cubic crystals the Cauchy relations can be 
written as an equality X = p between two Lame coefficients. The 
Cauchy relations apply if interactions between atoms (or colloids) 
are pairwise, central, sufficiently short ranged, and if thermal fluc¬ 
tuations are neglected. The role of thermal fluctuations was clar¬ 
ified by Squire et al . 11 who calculated the extra contributions 
(beyond those found by Born and Huang 10 ) in the expressions 
for the elastic constants at nonzero temperature. These contribu¬ 
tions can be seen as the “entropic” part of the elastic moduli, and 
they explain why hard-sphere crystals do not satisfy the Cauchy 
relations 12 . 

The paper is structured as follows: In section 2 we collect some 
results from known theories, in particular from the non-thermal 
dynamical (Hessian) matrix. We derive some properties of the 
dynamical matrix which allow comparison with the fully thermal 
case. The numerical data for the hard-sphere crystals is then pre¬ 
sented in section 3 in one, two and three dimensions. Finally, in 
section 4 we discuss the observed differences between the two 
cases. 

2 Effective interactions and the central 
force model 

The simplest (Cauchy) model of elasticity in the physics of solids 
supposes that interactions are pairwise and central between all 
particles in the solid, deriving from a scalar potential, see for ex¬ 
ample § 29 of the book by Born and Huang 10 . In this section 
we study the elastic properties of such systems, in order to better 
understand where the numerical data demonstrate the presence 
of non-central potentials. 

In this model, the total potential energy of a configuration is 
a function of all particle positions iq. It is expressed as a sum of 
central pair-potentials gy, 

*(M) = I«ij((ri-rj) 2 /2) (1) 

(«) 

The sum runs over all different pairs of particles. To simplify the 
algebra and avoid square roots in the calculations, the gy are func¬ 
tions of the distance squared. The functions should be identical 
for symmetry-related pairs of particles, but otherwise each func¬ 
tion is independent. 

In order to unify the notation with the numerical simulations 
of later sections, we assume already here that the particles are in 
a crystalline configuration (hexagonal or face-centered cubic, de¬ 
pending on the number of spatial dimensions). Throughout the 


paper, i,j,m,n denote tuples of integer indices on the Bravais lat¬ 
tice and serve to identify the sites on the lattice. Spatial indices 
are denoted by Greek letters. The spatial (Euclidean) positions of 
lattice sites are denoted by vectors Ri. They are multiples of the 
lattice spacing do. The numerical, and optionally also the theoret¬ 
ical calculations, are periodic in space. All differences between 
lattice indices, i —j, or lattice positions, Ry := Rj_j = R[ — Rj are 
thus understood modulo the periodicity, the result being mapped 
back into a centered copy of the periodic box. 

2.1 Description in terms of local displacements 

We now calculate the Hessian matrix of the potential en¬ 
ergy <5>({r}) around the crystalline reference state. This amounts 
to small-amplitude displacements, which are non-thermal in their 
nature. The reference state has inversion symmetry (hexagonal in 
2D, face-centered cubic in 3D), which imposes that the partial first 
derivatives vanish when evaluated at the reference positions: 

£<«) = o ® 

or \a 

for all i and a. From Eq. (2) it follows that the increase in poten¬ 
tial energy of the system due to local displacements is given to 
quadratic order by 

A<J> := <E>({r}) — <E>({R}) = ^ £ H a p(j-i)u ia u ip (3) 

a,1 3,i,j 

where Ui = iq — R[ is the displacement vector of particle i. The 
Hessian matrix consists of second derivatives, evaluated at the 
reference positions, 

(4 > 

Start with the first derivative 

If (W) = Ete-*i)*« ® 

' ifi 

and obtain H aji (j — i) for i # j: 

= -(R ia -Ria)(Rip-Rjp)gij-S a pg^ ( 6 ) 

We abbreviate our notation by introducing Rj_i = Rj — R[ and 
gj-i = gy and obtain 

H ap (n) = -R na R n p g" (R 2 J 2) - 8 ap g' n (R 2 J2). (7a) 

At j = i, the matrix is constrained by the fact that the sum over all 
n vanishes, 

HapW = - Y, H ap (n). (7b) 

n^O 

Let us return for a moment to the meaning of this Hessian ma¬ 
trix. As was mentioned in the introduction, it is also the dynam¬ 
ical matrix of particles (of unit mass), in the sense that their dis¬ 
placements oscillate around the reference state according to the 
equation 

“ia = -£#a/3(j-i) u j,8- (8) 

i,P 


2 
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In words, the interaction matrix H a p( i-j) governs the force ex¬ 
erted on particle i by particle j. For every separation vector n 
between two particles, H a p( n) is a dxd- matrix. We will find it 
important to analyse the forces and displacements parallel and 
orthogonal to the separation vector n. In particular in two dimen¬ 
sions, we can use the rotated orthonormal basis (R n ,o), where 
the hat denotes a normalised vector, to define the rotated compo¬ 
nents of the matrix H{ n). We find 



:= R^(n)R„ = ^"(R^/2) 

(9a) 

H±( n) 

:=d T H(n)d = -g , n (Ri/2) 

(9b) 

Ffsym (n) 

:= (R^(n)6 + 6 r ff(n)R„)/2 = 0 

(9c) 

^fasym( n ) 

:= (R^//(n)o —o r //(n)R n )/2 = 0 

(9d) 


The components #y and H_\_ control those components of the force 
which are parallel to the displacements, respectively both parallel 
to n (T/||) or both orthogonal to it (#J. The (a)symmetric parts 
govern those components of the force which are orthogonal to 
the displacements they result from. They are found to vanish 
in the present model, meaning that displacements can generate 
only forces which are parallel to them. This property appears as a 
result of the potentials being central, and it will be tested on the 
numerical data in the sections below. 

We can also examine the results in the non-rotated frame. In 
particular, for the off-diagonal terms we find 

Hxy(n) 4—X„y„g^(^/2). (10) 

which is an odd function with respect to X n and to Y n . Also this 
symmetry will be subject to comparison with the numerical data. 

2.2 Long-wavelength limit 

The elastic tensor for the model (1) has been calculated by Born 
and Huang 10 (§ 29). In our notation, their result for the elastic 
tensor reads 

Axct/It = X77 ^^n,a^n,cr^n,jS^n,T g n (^n/2) • (11) 

zv n 

Here, N is the number of particles, and V is the volume of the 
(periodic) box. We use an overbar to distinguish this tensor from 
the true (thermal) elastic tensor. Expression (11) is symmetric 
under all index permutations, which is higher than the symmetry 
required for the elastic tensor. This observation leads us directly 
to the Cauchy relation: The general form of the elastic tensor is 
given by group theory, for example in a two-dimensional hexago¬ 
nal lattice, 

Qx/3ctt = + m(<W^/3t + 3cj) j (12) 

with the Lame coefficients A and p. This expression becomes com¬ 
pletely symmetric under index permutations only if the Cauchy 
relation A = p is satisfied. The situation is similar in three- 
dimensional cubic lattices, only that there is another term in 
Eq. (12) with the anisotropy coefficient as prefactor. There, the 


Cauchy relation reduces three independent coefficients to two. In 
Section 3 we will measure the Lame coefficients from numerical 
data and test directly whether the Cauchy relation is satisfied. 

When we determine the Lame coefficients from Eq. (11), we 
find indeed that they are identical, 

= 03) 

In the last step we made use of Eqs. (9a) and (9b). 

We must mention that the calculation by Born and Huang 10 is 
done with a stress-free reference state in mind. In model (1) we 
do not know a priori whether the stress vanishes in the reference 
state or not, and for the hard-sphere crystals in the sections be¬ 
low, we know that they stay crystalline only if they are sufficiently 
compressed by the periodic box. We therefore need to consider 
the general case of a non-vanishing stress in the reference state. 
The effect of this stress is that the elastic tensor C, which char¬ 
acterises the quadratic increase of free energy with strain, is not 
identical to the elastic tensor A which characterises the motion of 
long-wavelength vibrations and thermal excitations. 13 The two 
tensors are related by 

^a/3cJT * = Qx/3<tt + (14) 

where T denotes the stress tensor in the reference state. To obtain 
the stress tensor, we start with the force between two particles, 

fy=-(R i -R j )^ j . (15) 

The virial part of the stress tensor is thus 

f «=-Ai( R i- R jUj,T 

(ij) 

= ov ^Rn,oRn,xg n - ( 16 ) 

n 

This expression has been given also by Born and Huang 10 . 

It is desirable to derive expression (11) directly from a long- 
wavelength limit of the interaction matrix H. This calculation 
cannot be done at this point, since we require thermal fluctuations 
which then yield the tensor A. From there we can conclude on 
the tensor C. We therefore defer this calculation to section 3.2.2 
below, where we will do it on the matrix D. When we replace 
D(n) by H(n)/kT in the result (30) there, we find the equivalent 
of the elastic tensor A, 

:= ~ ^7 ( n ) 

= 2 y T (^n,a^ n ,/3 &n + ^a/3 £n) 5 (17) 

which is entirely compatible with Eqs. (11), (14) and (16). 
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2.3 Description in terms of global deformations 


In the above section we expressed the elastic constants in terms 
of a long-wavelength limit of the matrix H. This matrix, as it 
involves second derivatives of the potential energy only is non- 
thermal in nature. At least for the elastic constants, it is possible 
to include thermal effects 11 . We introduce a parameter for global 
deformations of the system, the strain q and perform a Taylor 
expansion of the Free Energy F{rf)m —kT\n(Z(r\)) about zero 
strain (q = 0). As the whole box is changed parametrically, we 
could speak of a “zero-wavevector” (|Q| =0) method instead of 
the long-wavelength limit |Q| —>► 0. The resulting stress tensor and 
elastic tensor are expressed by the usual Gibbs-weighted average 
over configurations, 11 


_ 1 dF 

aP ~Vd^ 


(0) 


= -k T ^S a p + U E gy (r{ j/2) r ij a r ij fj ) (18) 

V V \(ij) 

1 d 2 F N 

Ca ^ X = Vd^d^ {0) = 2kT y 5 ^ 5 fi° 

v K1 (y) 

+ y ( L #ij ( r jj/ 2 ) r ij,a r ij,/3 r ij,ci r ij(19) 

v x (ij) 


The function Ccc (a,b) := (ab) - (a)(b) is a cumulant-like cross¬ 
correlation function. The terms on the right-hand side of Eq. (18) 
are called the kinetic and the virial terms. The terms in Eq. (19) 
are called the kinetic term, the fluctuations term, and the Born 
term. 11 


3 Numerical results from Hard sphere sim¬ 
ulations 

We performed high statistics event-driven molecular dynamics 
simulations of N hard spheres in one, two and three dimensions 
{d = 1,2,3). In one dimension the particles were confined to a 
line. In two dimensions we simulated a hexagonal crystal within 
a box adapted to the Bravais lattice. In three dimensions we per¬ 
formed simulations of face-centred cubic crystals, again within a 
box adapted to the Bravais lattice. In all dimensions the simula¬ 
tion box was continued periodically. During a simulation the total 
momentum and the (kinetic) energy were conserved. The densi¬ 
ties were always so high that no particle interchange occurred 
over the simulation time, so that particle diffusion and defects 
can be neglected in the data analysis. The Bravais lattice thus 
defines the reference positions of the spheres, and these coincide 
with the average positions. 

For all numerical data the units are chosen such that the par¬ 
ticles have unit diameter, unit mass, and that kT = 1 with T the 
temperature. 

During the simulation we took a series of K data recordings. 
For each recording, we measure the deviations Ui := ri — Rj of 
the instantaneous particle positions from their reference positions 
and calculate the dNxdN correlation matrix 

1 K 

G(i,a)(j,j3) : =^E u ia u jf) (20) 

A k= 1 

The matrix has a simple physical interpretation from linear re¬ 
sponse theory: It determines the vectorial displacement at j due 
to a force at i. We will thus call it a Green function in the follow¬ 
ing. In order to reduce statistical noise, we average elements of G 
which are related by the translational symmetry of the lattice, 

G ap(. n ) '-= ( 21 ) 


The virial term in (18) and the Born term in (19) resemble 
much the expressions we found for stress (16) and elastic ten¬ 
sor (11) from model (1). Only, there we evaluated at the refer¬ 
ence positions, and here we average over positions around these 
reference positions. Of course, if the position distribution is very 
narrow, the values obtained from both can be close. Both kinetic 
terms and the fluctuation term are missing, however. The absence 
of these terms is where the tensor C has its higher symmetry from, 
and why we obtained the Cauchy relation for it. 

The elastic tensor (19) includes all thermal effects, because we 
take the second derivative of the free energy. A similar derivative 
of the potential energy (1) would simply yield again the Born-like 
term in (11). It would be very nice to add thermal effects to the 
interaction matrix, just by replacing the potential energy by the 
free energy in the second derivative with respect to the reference 
positions Ri. Unfortunately, the reference positions are not pa¬ 
rameters, they appear by a spontaneous symmetry breaking, and 
we therefore do not know how to do such a derivative. For the 
moment, we can only work with the imperfect expressions of the 
above sections and compare them with full numerical data. 


The sum runs over all sites of the Bravais lattice. We could also 
average over the rotational and inversion symmetries, but we pre¬ 
fer to see the symmetries appear from the data in order get a idea 
of the statistical errors; the translational average appears to be 
sufficient for noise reduction. A visual representation of the ele¬ 
ments of G(n) in two dimensions is given in Fig. 1. It is a rather 
structureless object, and it depends on the size and shape of the 
periodic box due to the logarithmic nature of two-dimensional 
Green functions. In the figure we plot the elements of G(n) in a 
rotated coordinate system with one basis vector aligned with the 
vector R n . We note in particular the interesting physical structure 
that occurs in the off-diagonal elements. These elements of the 
response function vanish along high symmetry directions: a force 
along these symmetry directions gives a purely parallel displace¬ 
ment. This gives rise to a set of radial white lines in the lowest 
panel of Fig. 1. 

We then use matrix algebra to numerically invert the large ma¬ 
trix G to produce the effective interaction matrix D. More pre¬ 
cisely, we use the translation-averaged dNxdN matrix 

G ap(^i) : = G ce/3(j — »)• C 22 ) 


4| 
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Fig. 2 The effective interactions D(j-i) in a one-dimensional chain of 
N = 100 impenetrable rods. Beyond nearest neighbours the interactions 
are zero to within statistical noise. Larger blue symbols are positive 
values, smaller red are symbols negative values. K = N x 10 7 recordings, 
obtained from three independently simulated trajectories. 
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Fig. 1 The two-dimensional Green function in a periodic 32x32 
hexagonal lattice. Plotted are the components of the 2x2 matrices 
G aj8 (n) in a rotated orthonormal basis (R n ,6): Gy = R^GRn, G± = 6 r Go, 
G sym = (RnG6 + o r GR n )/2. The antisymmetric off-diagonal terms vanish 
to within statistical noise (±10~ 9 ). In the lowest panel the white lines 
correspond to a mechanical response parallel to an imposed force. 


and take its Moore-Penrose pseudoinverse because we must avoid 
inversion of the zero eigenvalues in G which come from momen¬ 
tum conservation. The inverse is denoted by D a p( i,j). In the 
following, we will analyse the properties of the resulting 2x2 ma¬ 
trices 

D a p(n):=^D a p(i,i + n ). (23) 

3.1 One-dimensional system 

Previous work has led to a detailed understanding of the one¬ 
dimensional problem. 14 We do not treat it in detail, but rather use 
it to check some of the chain of data analysis. The result is very 
simple - that the effective interaction in one dimensional fluids is 
limited to nearest neighbours; for larger separations D is zero. We 
used our code for simulation and data analysis and confirmed this 
results. Beyond the nearest neighbour, the effective interaction 
D(j — i ) in Fig. 2 falls to zero within statistical noise. Note that 
the relative statistical noise in this plot is at the level of 10 -5 , 
showing the high precision of the simulations. Given the excellent 
results found in one dimension, we can feel confident that the 
very different results found in two and in three dimensions are a 
result of differing physics, and not problems related to systematic 
or statistical errors in the data sets. 

In the data of Fig. 2 we find the ratio between the two nonzero 
values to be D(0)/D(1) = -2.0± 10 -5 . This is precisely the ratio 
that is expected from a discretisation of the Laplace operator in 
terms of nearest neighbours. We thus found the expected result, 
namely that the inverse of the (static) Green function is the corre¬ 
sponding differential operator of one-dimensional (static) elastic¬ 
ity. 

3.2 Two-dimensional crystal 

We performed two-dimensional simulations in systems with N = 
LxL hard disks, L varying from 10 to 100 with fixed surface frac¬ 
tion (j) = 0.85. We expect that the elements of the effective dynam¬ 
ical matrix D come from the local physics (and available entropy) 
occurring near each particle and are not the results of a non-trivial 
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Fig. 3 Effective interactions as a function of the inverse system size l/L. 
Different plot styles correspond to different neighbour layers, starting 
from 0 at the top, down to 7. Fewer layers are plotted for L = 10. 



Fig. 4 Effective interactions in a hexagonal hard-disk crystal with L = 32. 
Larger blue symbols are positive values, smaller red symbols are 
negative values. K = 2Nx 3.5xlO 6 recordings, obtained from 
336 independent trajectories. Surface fraction </> = 0.85. These data 
required around 300 000 CPU core hours. 


propagation of boundary effects down to the microscopic scale. If 
the physics is local, we then expect that the values of the elements 
of D vary very little with changes in the system size. To test this, 
we plot in Fig. 3 the evolution of tr(D(n)) with the system size L 
and find that the lines for small separations n are remarkably sta¬ 
ble. This is true for all matrix elements of D a p, not only for the 
trace. Since we are performing simulations within a periodic box, 
it is clear that we should only be evaluating elements for sepa¬ 
rations such that periodic copies do not contaminate the result. 
Thus we will only quantitatively analyse D for separations which 
are smaller than L/2. 

The correlation function G, which is measured in the simula¬ 
tions and displayed in Fig. 1, is subject to strong finite size cor¬ 
rections due to the logarithmic nature of two-dimensional Green 
functions. The components of the matrix D, however, do not 
evolve for small values of l/L. From the curves in Fig. 3 we 
decided to use a fixed value of L = 32 for further simulation in 
order to collect the highest possible statistics, while being able to 
resolve effective interactions out to a large distance. 

The resulting matrix elements are plotted in Fig. 4. As we did 


already in Fig. 1, we rotated the 2x2 matrices D a p( n) for each 
separation vector R n . We used the same rotated orthonormal ba¬ 
sis and the same notation as we used in Eqs. (9). The interaction 
matrix depicted in Fig. 4 does not vanish beyond the first layer of 
neighbours - very differently from its unidimensional counterpart 
in Fig. 2. We here observe that D(n) is not sparse. We used very 
high statistics in this plot to be sure that the interaction data is 
well separated from statistical noise - We find it to be the case for 
the first six layers of neighbours. The noise level can be read off 
from the values of D asyin which are zero in noiseless data, finding 
values below 2x 10 -2 . Beyond the 10th layer of neighbours, signal 
and noise cannot be distinguished anymore. Another measure of 
the statistical noise is the spread of those symbols which are equal 
by symmetry of the lattice (smaller than the linewidth in the plot). 
Both measures for the noise were observed to decrease while we 
accumulated more and more data during the simulations. 

We tried to characterise the decay of the effective interactions 
with particle separation, which we can resolve out to the eighth 
layer of neighbours. We tried fitting the data with both exponen¬ 
tial and power law decays, plotting both particle separation and 
layer number for the abscissa. No fit seemed totally convincing 
with our data. If one insists on fitting with a power law ||R n || a , 
one finds a ^ — 6 ±0.5. At the moment we do not have analytic 
arguments that would predict the functional form of the decay. 

3.2.1 Two-dimensional visualisation 

We now give some more details of the data presented in Fig. 4. 
In particular, we like to visualise those points which are charac¬ 
terised by being within same neighbour level, but which differ 
with respect to the hexagonal symmetry. 

We can study the 2x2 matrices D(n) in different frames of refer¬ 
ence. The first frame of reference is the orthonormal basis (R n ,o) 
already used above. In Fig. 5 we plot the parallel, the orthog¬ 
onal, and the off-diagonal component of the matrix. All three 
panels are invariant under rotations of 60 degrees. Under mirror¬ 
ing, however, only the first two are invariant while D s ym changes 
sign. We did not plot D aS ym because it contains only noise. 

A second manner of examining the data is in the fixed Carte¬ 
sian frame, which is the same for all n. In Fig. 6 we study the 
trace of the matrix, tr(D), the first component D xx of which we 
subtracted half the trace, and the off-diagonal element Dxy , sym¬ 
metrised. The component D yy — tr(D)/2 is not plotted as it is sim¬ 
ply the negative of D xx — tr(D)/2. Only tr(D) still has hexagonal 
symmetry. The two bottom panels show (skew) symmetry with 
respect to mirroring about the v-axis and the y-axis. We found the 
same skew symmetry in the Hessian matrix, Eq. (10). The special 
choice of component combinations in Fig. 6 is motivated by the 
continuum limit which we discuss now. 

3.2.2 Hexagonal continuum limit 

The matrix D a p (n) is the inverse of the Green function of static 
elasticity and is thus expected to be related to the correspond¬ 
ing differential operator - or rather its hexagonal discretisation. 
Continuum elastic theory of pre-stressed materials gives the dif¬ 
ferential operator (Sec. 3 of Ref. 13) 

&ao = — (24) 
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Fig. 5 Matrix elements of D(n) in rotated frames of reference, for Fig. 6 Some linear combinations of matrix elements D a p( n) for 

several R n . The center of the hexagon corresponds to R n = 0. The several R n . The antisymmetric off-diagonal (D*y —D yx )/2 vanishes to 

antisymmetric off-diagonal D asym vanishes to within statistical noise. within statistical noise. 

The reflection anti-symmetry of D sym imposes that the elements of D sym 
are zero for the first two neighbour layers. All elements of D sym are small 
compared to Dm and D ± . 
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where we have to remember that the two elastic tensors A and C 
coincide only in systems whose reference state is stress-free. The 
hard-sphere crystal in our numerical simulation is compressed by 
the periodic boundary conditions, leading to the isotropic stress 
Tp z = —P5p T with the pressure P. Consequently the two tensors 
A and C are different. Assembling Eqs. (12), (14), we find the 
differential operator to be 


^^(jU-P)5 a/3 V 2 -( k + fl)d a d p . (25) 


For these second derivatives, we can produce simple discretisa¬ 
tions in terms of finite differences on nearest neighbours. In par¬ 
ticular, the combinations of matrix elements used in Fig. 6 are 
simple combinations of second derivatives,* 
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(26b) 


^cy + ^yx ^ _dxdy + dyd, = (26c) 

2 2 

We see that these stencils are qualitatively reproduced in the cen¬ 
tral parts of the panels in Fig. 6. However, Fig. 6 shows nonzero 
interactions also beyond the first layer of neighbours. Another dif¬ 
ference is that for tr(D) we find a middle value which is —5.6 times 
the values found on the first neighbours (instead of —6). When 
comparing the values of the first neighbours among themselves, 
the stencils are again well reproduced, with variations as small 
as ~ 10 -5 in all three panels of Fig. 6. The skew symmetry of the 
operator d x d y is found again in all values of Dxy(n). This symme¬ 
try imposes the horizontal and the vertical white lines (zeros) in 
the third panel of Fig. 6. 


The matrix D a p (n) encodes the full dispersion curves that are 
required to determine the spectral properties of the fluctuations. 
In particular, one can extract the elastic tensor A from the long- 
wavelength limit. The usual way to do this calculation is to 
perform a discrete periodic (fast) Fourier transform on the ma¬ 
trix G a p( i,j) and to observe that the translation invariance ren¬ 
ders the result block-diagonal. For each reciprocal vector Q m we 


* These stencils are obtained by writing the operator in question as V r MV, with a 
corresponding matrix M (for example, M t j = 8^ for V 2 ). The matrix is then decom¬ 
posed into a suitable sum of terms e/ej, where the e,- are the six hexagonal unit 
vectors obtained from rotation by 60°. The result is thus composed of terms such as 
ef V, which are directed derivatives between the center and a first neighbour. Each 
of these terms is thus immediately translated into a finite difference between them, 
resulting in the weights noted in the stencils of Eqs. (26). 
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Fig. 7 Dispersion curves used to obtain the elastic moduli in the 
limit |Q| 0, here for L = 100. The dotted lines indicate values from a 

different numerical method working at |Q| = 0. 


have one Fourier transformed dxd matrix 

G ap (m):=N'Ee- iQm -*'G a i i ( n) (27) 

n 


= (u a (m)up (m) 


(28) 


The overline denotes a complex conjugate. To leading order, this 
matrix scales as |Q m | -2 , such that we can extract the elastic tensor 
as the long-wavelength prefactor to this scaling, 


AaGpTQoQz — 


kTN 2 
V 


lim 

|QmK0 


[ 6 ( m ) ‘W 

|Qm | 2 


(29) 


In practice, rather the eigenvalues of the inverse matrix are plot¬ 
ted, see Fig. 7. The long-wavelength limit is then done by extrap¬ 
olation by eye. The upper curve corresponds to the longitudinal 
waves, converging to A XXX x = (M — P) + (A + /i) in the limit |Q| -> 0. 
The lower curve gives the transverse value A*y*y = /i — P. 

The same limit can be done using the matrix D instead of G. 
The block structure helps to invert the large matrix, the inverse 
is again block-diagonal. The dxd blocks are simply the Fourier 
transforms of the matrices D a p( n): [G{m)~ l ] a p = D a p(m)/N 2 . 
For the long-wavelength limit we expand the exponential in the 
definition of the Fourier transformation (27) around Q m = 0. The 
sum over the constant term vanishes, £ n D a ^(n) = 0. The linear 
term vanishes due to symmetry, such that we obtain the quadratic 
term as the leading one. This conveniently matches with what we 
require for the limit in Eq. (29). Factorising the unit reciprocal 
vectors, we obtain for the elastic tensor 

Aao(3z = ~ 2 y ( n )- (30) 


To determine the Lame coefficients, we can either choose individ¬ 
ual components such as A xxxx above, or we can do a “hexagonal 
average”, which amounts to reducing the rank-4 tensor to scalars, 
for example with 8 a p8 oz and with 8 aa 8p z . This gives a 2x2 sys¬ 
tem of equations for (fi -P) and (A + fi), which when solved be- 
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comes 



Fig. 8 The summands of Eqs. (31) and (32). 



Fig. 9 Convergence of the sums in Eqs. (31) and (32). 


> x - p = -l^Y-'L R n( 3z M n )+£||(n)) (31) 

A+ai = n)-D||(n)) (32) 

The summands are displayed in Fig. 8 . They decrease with in¬ 
creasing distance, as expected, and their sums converge. We visu¬ 
alise the convergence in Fig. 9 where we added up the summands 
layer by layer. The final values ar e fi-P = 373 and A + /x = 553. 

In order to extract the Lame coefficients, we need also the value 
of the pressure P. During simulations, we tracked the average flux 
of linear momentum, which is a mechanical definition of the pres¬ 
sure tensor. We further measured the response of this tensor to 
small deformations of the simulation box, which presents another 
(faster converging) way to calculate the elastic moduli. Notice 
that this method corresponds directly to “|Q| = 0”, no limit has to 
be taken. We obtain the values ^ 

P = 34.4, ii= 405, A = 143.5. (33) 

They are used in Fig. 7 to check the consistency of the different 
methods. Also the values extracted from Eqs. (31) and (32) are 
reasonably close. 

3.3 Three-dimensional crystal 

We also performed simulations in three dimensions, collecting 
high statistics data on a system of dimensions L= 15. Here, the 
result is similar to two dimensions, the effective interaction does 
not vanish beyond the layer of nearest neighbours but shows a 
rapid decrease. We adopt the same route of a rotated frame of ref¬ 
erence to plot in Fig. 10 the matrix elements of the 3x3 matrices 
D a p( n). We choose an orthonormal basis (R n , 6 i, 02 ), in which 
the matrix comprises the following blocks: 



As before, D y = RjDR n , but is a 2x2 matrix, and v, w are two- 
dimensional vectors. As the basis vectors 6 i, 62 can be chosen with 
an arbitrary rotation about the vector R n , we are only interested 
in invariants of D^,v,w under this rotation. For the matrix D j_ 
we thus plot in Fig. 10 the trace of Dj_, its anti-symmetrised off- 
diagonals D _\_ ?asym and some notion of its determinant. For v, w we 
plot the modulus after (anti)-symmetrising, D V:sy m := ||v + w||/2, 
^v,asym ||v — w||/2. The lowest panel in Fig. 10 allows one to 
estimate the noise level for the given statistics. 


t These values are compatible with those given in Ref. 12. We applied the box- 
deformation method also to the value 0 = 0.863714 which is given in that reference, 
and our implementation reproduces exactly the given values for the elastic moduli. 
Concerning the numerical values of elastic moduli, there was a disagreement in the 
literature between Ref. 15 and Ref. 12, see also subsequent publications. Given that 
we implemented both the fluctuation method and two deformation methods (in a 
second one the spheres are deformed instead of the box) which all give the same 
result, we think that we can resolve the debate in favour of Ref. 12. 
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Fig. 10 Effective interactions in a FCC crystal of N = 15 3 hard spheres. 
Larger blue symbols are positive values, smaller red symbols are 
negative values. K = 3iVx 6.8xl0 4 recordings, obtained from 
424 independent trajectories; volume fraction 0 = 0.57. These data 
required around 230 000 CPU core hours. 


If again one insists on an algebraic fit D(n) °c ||R/do|| a > one 
finds a « — 8 ± 1. In three dimension the data display less scatter 
when plotted in terms of neighbour layer, rather than separation. 
This is a non-trivial geometrical effect which is already visible in 
two dimensions: In the first panel of Fig. 5, the values on the 
diagonals of the hexagon are larger than the values at the same 
distance (even larger than those on the same neighbour level). 
This creates within every neighbour level (or distance) a tendency 
which is opposed to the general trend. 

4 Discussion and conclusions 

In section 3 we have presented a number of properties, some of 
which are compatible, others are incompatible with the expecta¬ 
tion elaborated from the central model in section 2. Let us now 
scrutinise these properties and see to what extend we can under¬ 
stand the differences. 

4.0.0.1 Elastic moduli: Let us start with the elastic moduli 
and the Cauchy relation. The values we found for the Lame co¬ 
efficients in Eqs. (33) show that the two-dimensional hexagonal 
crystal is far from being “Cauchy”, p being nearly three times as 
large as A, where the Cauchy relations would require them to be 
equal. This finding is neither new 12 nor surprising, as shows the 
calculation by Squire et al 11 , which we recapitulated in Sec. 2.3. 
One concludes that the non-thermal central force model is inad¬ 
equate for understanding the elastic moduli of the hard-sphere 
system, that the kinetic term and probably more important the 
fluctuations term are not small, and finally that we can indeed ex¬ 
pect the hard-sphere system to be dominated by thermal effects 
and to see important deviations also in the interaction matrix. 

4.0.0.2 Pressure: As a side remark, we find it an interest¬ 
ing question whether one can determine the pressure in a sys¬ 
tem by measuring displacements and their correlations. 16 Equa¬ 
tion (16) shows that this is indeed possible in the non-thermal 
model, where the pressure is minus a diagonal term of the stress 
in Eq. (16), or P = — ^ L n ^n&n (^n/2) • The values of are trans¬ 
lated back to H_\_ using Eq. (9b). Does the same work also in the 
thermal case, using instead of ifj_? When we use the numer¬ 
ical data for Dj_, shown in Fig. 5, then we find the value -96 for 
this effective pressure - which is clearly unacceptable, even the 
sign is wrong. In fact, a close look on Eqs. (31) and (32) reveals 
that the sum actually calculates P + (A —p)/2 and equals the pres¬ 
sure only to the extent that the Cauchy relation is satisfied. 

4.0.0.3 Locality of the interaction: Independent of all ques¬ 
tions of how the data compare to a model, we found that the 
correlation matrix strongly depends on the system size, whereas 
its inverse does not (Fig. 3). Our intuition, in which we called 
the correlations a “Green function” and its inverse a “differen¬ 
tial operator” heavily relies on a clear-cut separation of size de¬ 
pendence. We find that the differential operator CD) is indeed 
“local”, but the precise notion of this locality is different in one 
and in higher dimensions. In one dimension we find the inter¬ 
action matrix to be restricted to the two nearest neighbours, as 
expected from a Laplace operator. In two and three dimensions 
the results are much more interesting: Despite the true interac¬ 
tions taking place to nearest neighbours, we find longer-ranged 
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effective interactions from the analysis of displacement correla¬ 
tions. These interactions are small compared to the ones between 
nearest neighbours, but they are not exactly zero and can be dis¬ 
tinguished from noise. It remains open what functional form the 
decrease in Figs. 4 and 10 follows, and whether an algebraic form 
would still allow to reasonably call the operator “local.” 

4.0.0.4 Theoretical interpretation of the dynamic matrix as 
a self-energy: The presented results open the question of a bet¬ 
ter theoretical interpretation than the one in Sec. 2.1. In particu¬ 
lar one would like to predict how the effective interactions decay 
with distance. It is clear that in a field-theory formulation of a per¬ 
turbed Gaussian theory the effective interactions D correspond to 
the self-energy in a Dyson equation. However, without an explicit 
perturbation scheme for hard sphere systems it is difficult to turn 
this remark into an explicit calculation. 

4.0.0.5 Nonlinearities in the interaction potential: Getting 
back to the central-force model of Sec. 2.1, it does not favour 
local or nonlocal interactions in its general form. If we want more 
information, we need to specify the functions gy. One choice is 
to take the true hard-sphere interactions, the pair potential being 
either zero or infinity. Around the reference state, the potential is 
zero, together with its first and second derivative. Therefore, the 
Hessian matrix vanishes, and also the stress and elastic tensors 
consist only in the kinetic terms. Here, the model is clearly out of 
its applicability. 

We learn from this remark that the model equations do not 
necessarily serve to infer the true interaction between particles 
from given data. But, what lesson do we want to learn from the 
model? Our aim was to contrast thermal with respect to non- 
thermal effects. Unfortunately, in the hard-sphere simulations, 
together with thermal effects, we introduced strong nonlineari¬ 
ties in the pair potential. Now we have a hard time to separate 
their influence from the thermal effects. (At least in the elastic 
constants, the nonlinearities play no role, which is reassuring.) 

Another convenient choice for the functions gy is a network 
of linear springs which have zero length at rest. When such a 
network is stretched, the particles may be found on a hexagonal 
reference lattice. In this system the Hessian matrix is a constant, 
and all statistical averages are Gaussian integrals which can be 
solved analytically via Wick’s theorem. The effective interactions 
are indeed found to be given by the Hessian matrix, D = H/kT, 
which in turn is given by the connections of particles. If we choose 
to have identical springs between nearest neighbours only, then 
also the matrix D is restricted to nearest neighbours. The same 
applies to second-nearest neighbours, and so on. One might think 
that this linear system will help to differentiate between thermal 
and nonlinear effects. However, when one evaluates the averages 
in Eq. (19), which contains all thermal effects, one finds that the 
elastic tensor is zero.^ The zero-length spring network is thus not 
good enough to separate nonlinear from thermal effects. 


$ This does not mean that the sound waves in this system have vanishing sound veloc¬ 
ity Sound velocities as well as phonons are encoded in the tensor A rather than C, 
and only the latter vanishes. The tensor A has additional contributions which are 
proportional to the stress of the reference state, see Eq. (14). 


If we insist on central pair potentials, the zero-length spring 
network is the only linear one we know of. Of course, other ar¬ 
bitrary linear models can be created by choosing a potential of 
the form (3), i. e. exactly quadratic in the displacements - one 
can even take the measured matrix D as the matrix in this expres¬ 
sion, which amounts to the shadow system of Ref. 7. However, 
doing so one relaxes the constraint of the central potential and 
includes knowledge about the reference positions in the poten¬ 
tial. Moreover, by construction the shadow system has the same 
displacement correlations, the same interaction matrix, and the 
same elastic constants as the original system. Only higher-order 
correlations differ. We are thus back to the question what do we 
want to learn from comparing with a model? Should it contain 
information about the reference positions, thus about the sponta¬ 
neous symmetry breaking or not? 

Another possible approach to a model would consist in ex¬ 
tracting the “best” central-force model from the numerical data. 
Wanted is a discrete function g n of which the first and sec¬ 
ond derivatives are given by gj, = —Dj_(n) and g" = (Dj_(n) - 
D||(n))//?n> respectively - in the hexagonal case, in analogy to 
Eqs. (9a) and (9b). The first derivative can be read off Fig. 5, 
and the second one is given up to a factor in the lower panel 
of Fig. 8. We did not try to perform this extraction and do not 
know whether the problem has a unique solution or approxima¬ 
tion methods are required. 

4.0.0.6 Non-central forces: The main lesson we can learn 
from the comparison with the model is probably that the effec¬ 
tive interactions contained in D are non-central: It is striking that 
H sym (n) is identically zero for all n, whereas D sym ( n) has nonzero 
values. The thermal interactions are thus genuinely non-central. 
We note that the off-diagonal term D S ym(n) is smaller in ampli¬ 
tude than the terms on the diagonal, D y and D_\_, so that the non¬ 
central nature can be considered a correction to the dominant 
terms. 
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